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ABSTRACT 

Context. Models of the young solar nebula assume a hot initial disk such that most volatiles are in the gas phase. Water emission 
arising from within 50 AU radius has been detected around low-mass embedded young stellar objects. The question remains whether 
an actively accreting disk is warm enough to have gas-phase water up to 50 AU radius. No detailed studies have yet been performed 
on the extent of snowlines in an accreting disk embedded in a dense envelope (Stage 0). 

Aims. Quantify the location of gas-phase volatiles in the inner envelope and disk system for an actively accreting embedded disk. 
Methods. Two-dimensional physical and radiative transfer models have been used to calculate the temperature structure of embedded 
protostellar systems. The heating due to viscous accretion is added through the diffusion approximation. Gas and ice abundances of 
H 2 O, CO 2 , and CO are calculated using the density dependent thermal desorption formulation. 

Results. The midplane water snowline increases from 3 to ~ 55 AU for accretion rates through the disk onto the star between 1(L 9 - 
1(L 4 M q yr _1 . CO 2 can remain in the solid phase within the disk for M < 1(L 5 M Q yr~‘ down to ~ 20 AU. Most of the CO is in 
the gas phase within an actively accreting disk independent of disk properties and accretion rate. The predicted optically thin water 
isotopolog emission is consistent with the detected Hj s O emission toward the Stage 0 embedded young stellar objects, originating 
from both the disk and the warm inner envelope (hot core). An accreting embedded disk can only account for water emission arising 
from R < 50 AU, however, and the extent rapidly decreases for M < 10~ 5 M e yr _1 . Thus, the radial extent of the emission can be 
measured with future ALMA observations and compared to this 50 AU limit. 

Conclusions. Volatiles such as H 2 O, CO 2 , CO, and the associated complex organics sublimate out to 50 AU in the midplane of young 
disks and, thus, can reset the chemical content inherited front the envelope in periods of high accretion rates (> 10 -5 M e yr~'). A hot 
young solar nebula out to 30 AU can only have occurred during the deeply embedded Stage 0, not during the T Tauri phase of our 
early solar system. 
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1. Introduction 


The snowlines of various volatiles (sublimation temperature 
E su b ^ 160 K) play a major role for planet formation. Beyond 
the snowline, the high abundances of solids allow for efficient 
sticking to form larger bodies, which is further enhanced by the 
presence of ices (e.g. Stevenson & Lunine 1988; Ros & Johansen 
2013). Extensive studies have investigated the snowline in pro¬ 
toplanetary disks around pre-main-sequence stars similar to the 
nebula out of which supposedly the solar system formed (e.g. 
Lissauer 1987; Pollack et al. 1996). In such models, the water 
snowline is located at a few AU radius. It is thought that the early 
pre-solar nebula was hot (> 1500 K) such that both volatiles and 
refractories (7) uh > 1400 K) are in the gas phase out to larger 
distances (Cassen 2001; Scott 2007; Davis et al. 2014; Marboeuf 
et al. 2014). The evidence of such a hot solar nebula comes from 
the history of the refractories, however the volatile content of 
comets seems to indicate that a part of the disk remains cold 
(Bockelee-Morvan et al. 2000; Mumma & Charnley 2011; Pon- 
toppidan et al. 2014). The evolution of the snowline due to disk 
and star evolution and its accretion rate clearly affects the chem¬ 
ical composition in the region relevant to planet formation (e.g., 
Lodders 2004; Davis 2005; Oberg et al. 201 lb). The most rele¬ 


vant volatiles are the known major ice species: H 2 O, CCL, and 
CO. Observations (Meijerink et al. 2009; Zhang et al. 2013) and 
models (e.g., D’Alessio et al. 1998; Dullemond et al. 2007) of 
protoplanetary disks around pre-main sequence T-Tauri stars in¬ 
dicate that such disks are not warm enough to have gas-phase 
volatiles in the midplane beyond 30 AU as claimed in some early 
solar nebula models, and, for the case of H 2 O, a snowline of only 
a few AU is typically found. Higher temperatures at large radii 
could potentially be achieved, however, during the deeply em¬ 
bedded phase of star formation when the accretion rate is high. 
The question remains how hot can an embedded accreting disk 
be when the accretion rate is high (> 10 6 M 0 yr 1 , see Dunham 
et al. 2014 for a recent review). 

Significant progress has been made in identifying snowlines 
in protoplanetary disks in the later stages when the envelope 
has dissipated and on their location with respect to gas giant 
formation sites (e.g., Kennedy & Kenyon 2008; Pontoppidan 
et al. 2014). Direct observational evidence of snowlines of the 
major ice species toward protoplanetary disks around pre-main 
sequence stars rely on the chemical changes that occur when a 
molecule is absent in the gas phase. The most readily observed 
snowline is that of CO as inferred through spatially resolved ob- 
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servations of N 2 H + , whose gas-phase abundance is enhanced 
when CO is frozen-out (e.g. Qi et al. 2013). DCO + is also a 
tracer of cold gas at temperatures close to that of the CO snow¬ 
line (van Dishoeck et al. 2003; Guilloteau et al. 2006; Qi et al. 
2008; Mathews et al. 2013). Both tracers indicate CO snowlines 
at >30 AU for T Tauri disks to >100 AU for disks around Her- 
big stars. The water snowline has been inferred to be within a 
few AU from direct observation and modelling of mid-infrared 
water lines (Meijerink et al. 2009; Zhang et al. 2013). The CO 
snowline location with respect to those of water and CO 2 has a 
direct impact on the amount of water present in giant planets at¬ 
mospheres (Oberg et al. 201 lb; Madhusudhan et al. 2011; Moses 
et al. 2013). 

Snowlines in the early embedded stages of star formation can 
be much further out, however, since the stellar accretion pro¬ 
cess through the disk onto the star begins at the time that the 
disk itself is forming. The gaseous volatile reservoir in embed¬ 
ded disks is affected by their formation process, which results in 
emerging protoplanetary disks having different chemical struc¬ 
tures depending on initial cloud core parameters (Visser et al. 
2011). A related question centers on whether these volatiles are 
‘inherited’ or ‘reset’ during the planet formation process (Pon- 
toppidan et al. 2014). The ‘reset’ scenario refers to the chemical 
processing of ices as the gas and dust are exposed to elevated 
temperatures (> 40 K) during their voyage from the envelope 
to the disk. These temperature regions define the regions where 
both CO 2 and CO sublimate from the ice in the gas phase. Both 
species are however difficult to trace in embedded disks: CO 2 
because it lacks a dipole moment and CO because of confusion 
with the surrounding envelope. 

The effect of accretion may be most readily seen through the 
changes in the water snowline. Water is the major constituent 
of ices on the grains that facilitate planet formation (Gibb et al. 
2004; Oberg et al. 2011a) and a major coolant (Karska et al. 

2013) , and thus a key volatile in star- and planet formation. Most 
of the water is thought to be formed during the pre-stellar stage, 
and then transported through the envelope into the disk and plan¬ 
ets (Visser et al. 2009; Cleeves et al. 2014; van Dishoeck et al. 

2014) . The crucial step that is yet relatively unexplored is the 
processing of water during disk formation. 

The water vapor content around protostars is investigated by 
the ‘Water In Star-forming regions with Herschel’ key program 
(WISH, van Dishoeck et al. 2011). Due to the large beam of 
Herschel, a significant fraction of the detected water emission 
is from the large-scale envelope and the bipolar outflow (Kris- 
tensen et al. 2012; Herczeg et al. 2012; Mottram et al. 2014). 
More importantly, the outflowing water will escape the system 
and will not be retained by the disk. In order to determine the 
amount of water vapor associated with the inner envelope or 
disk, an isotopolog of warm water (H' s O) needs to be observed. 
Visser et al. (2013) reported a detection of the H[ x O line (1096 
GHz E u = 249 K)) with HIFI (de Graauw et al. 20l0) toward the 
embedded protostar NGC 1333 IRAS2A. However, they found 
that the line is still optically thick with an emitting region of 
~ 100 AU. Thus, it remains difficult to constrain the amount of 
water vapor in embedded disks through single-dish observations. 

Spatially resolved warm H' s O emission (excitation temper¬ 
ature T ex ~ 120 K Persson et al. 2014) has recently been de¬ 
tected from within the inner 50 AU radius of several deeply 
embedded (Class 0) low-mass protostars (e.g., Jprgensen & van 
Dishoeck 2010; Persson et al. 2012, 2013) and toward one high 
mass disk (van der Tak et al. 2006; Wang et al. 2012). The low- 
mass sources are very young objects whose envelope mass is 
substantially higher than the mass at small scales (R ^ 100 AU) 


(also denoted as Stage 0, Robitaille et al. 2006). The emission 
is expected to arise from 7)| ust > 100 K regions where ice sub¬ 
limates and the gas-phase water abundance is at its maximum 
(Fraser et al. 2001; Aikawa et al. 2008; Mottram et al. 2013). 
This inner region is, however, also where the disk forms (Larson 
2003; Williams & Cieza 2011; Li et al. 2014). Rotationally sup¬ 
ported disks have been detected recently around a few low-mass 
protostars (e.g., Tobin et al. 2012; Murillo et al. 2013; Ohashi 
et al. 2014). From the kinematical information, both Jprgensen 
& van Dishoeck (2010) and Persson et al. (2012) found that the 
water emission does not show Keplerian motion and concluded 
that it must be emitted from a flattened disk-like structure that is 
still dominated by the radial infalling motions. 

This paper investigates snowlines of volatiles within an ac¬ 
creting disk embedded in a massive envelope. The spatial extent 
and the water vapor emission are compared with the observed 
values toward three deeply embedded low-mass protostars. The 
thermal structure of an actively accreting disk is computed in¬ 
cluding the additional heating due to the energy released from 
the viscous dissipation. Most previous studies of the thermal 
structure of an accreting disk focused on the later evolution¬ 
ary stage of disk evolution where the envelope has largely dissi¬ 
pated away. Furthermore, they focused on the midplane temper¬ 
ature structure (e.g., Sasselov & Lecar 2000; Lecar et al. 2006; 
Kennedy & Kenyon 2008). The additional heating in the em¬ 
bedded phase shifts the snowlines of volatiles outward to larger 
radii than in disks around pre-main sequence stars (Davis 2005; 
Garaud & Lin 2007; Min et al. 2011). The details of the phys¬ 
ical and chemical structure of the embedded disk are presented 
in Section 2. Section 3 presents the snowlines location as func¬ 
tion of disk and stellar properties. The results are compared with 
observations and their implications on the young solar nebula is 
discussed in Section 4. Section 5 summarizes the main results 
and conclusions. 

2. Physical and chemical structures 


Table 1. Parameters for the embedded disk + envelope models. The 
varied parameters and the canonical values are indicated in boldface. 


Variable [unit] 

Description 

Value(s) 

rout [AU] 

Outer radius 

lO 4 

n„ [AU] 

Inner radius 

0.1 

Rcen [AU] 

Centrifugal radius 

50, 200 

M env [ M 0 ] 

Envelope mass 

1.0 

Mdisk [M 0 ] 

Disk mass 

0.05, 0.1, 0.2, 0.5 

Rdisk [AU] 

Disk radii 

50, 100, 200 

H 0 [AU] 

Scale height at 1 AU 

0.2 

r* [K] 

Stellar temperature 

4000 

M* [M 0 ] 

Stellar mass 

0.5 

L* [As] 

Stellar luminosity 

1,5, 15 

M [M 0 yr _1 ] 

Accretion rate 

J Q-4,-4.2,-4.55 - 6 - 1-9 


2.1. Physical structure 

A parametrized embedded disk (disk + flattened envelope) 
model is used to construct the density structure following Crapsi 
et al. (2008). The main parameters are disk mass (Mdisk) and disk 
radius (Rdisk)- A number of parameters defining the envelope and 
the disk are fixed and summarized in Table 1. The envelope mass 
is fixed at 1 M 0 , which is appropriate for the objects from which 
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the water emission has been detected (J0rgensen et al. 2009; 
Kristensen et al. 2012, and see Section 4.3). The mass distri¬ 
bution within a flattened envelope is more crucial for the tem¬ 
perature structure than the total mass. 

For the large-scale envelope, a flattened envelope due to ro¬ 
tation as described by Ulrich (1976) is adopted whose densities 
are given by the following equation: 



where /./ = cos 6, R cen the centrifugal radius, and r the spherical 
radius. The centrifugal radius defines the region in which the ma¬ 
terial no longer flows radially and enters the disk. Two centrifu¬ 
gal radii of R cen = 50 AU and 200 AU are explored. The low R cen 
corresponds to the small Keplerian disk toward L1527 (Ohashi 
et al. 2014) while the high R cen corresponds to corresponds to 
the maximum disk radius (~ 180 AU) observed toward a Class 
0 embedded low-mass YSO (Murillo et al. 2013). The two cases 
explore the effect of mass concentration in the inner envelope on 
water emission. For a given centrifugal radius, a particle follows 
a parabolic motion given by 


r l-fi/fip 
Rcen 1 — /Tq 


( 2 ) 


where po satisfies the condition above at every r and p. The outer 
radius of the envelope is fixed at r out = 10 4 AU with an inner 
radius of 0.1 AU where the dust typically sublimates (assuming 
a dust sublimation temperature, 7/ uh , between 1500-2000 K). 

An outflow cavity is then carved out from the envelope den¬ 
sity structure at > 0.95. Following Crapsi et al. (2008), the 
density inside the cavity is equal to that of the densities at r out . 
This creates a conical outflow with an aperture of 30° at large 
radii (semi-aperture of 15°). For a 1 M 0 envelope, gas densities 
of ~ 10 4 cirT 3 fill the envelope cavity, which is consistent with 
those observed toward YSOs (e.g., Bachiller & Tafalla 1999; 
Whitney et al. 2003). 

A flared accretion disk is added to the envelope density struc¬ 
ture. The density within the disk follows the power law depen¬ 
dence radially and has a Gaussian distribution vertically in z. 
as expected from a hydrostatic disk. The flared disk densities 
(Shakura & Sunyaev 1973; Pringle 1981; Hartmann et al. 1998; 
Williams & Cieza 2011) are described by 


Pdisk (R, z ) 


£o x (R/Rdisk) 1 
yf2^H(R) 


exp 


1 

2 



(3) 


where the scale height H is fixed to 0.2 AU (Hq) at 1 AU (Rq), 
7?disk the disk radius, and R the cylindrical radius. The radial 
dependence of the scale height is H(R) = R Hq/Rq (R/Rq) 2/1 
(Chiang & Goldreich 1997). Finally, the densities are scaled by 
a constant factor Zq such that the total disk mass within TGisk 
(R < 7?disk) is equal to the values in Table 1 and distributed 
within /?disk- In the case of R cen = 50 AU flattened envelopes, 
only Pdisk = 50 AU models are considered. The total gas density 
in the model is p — pdisk + p e nv with a gas-to-dust mass ratio of 
100 . 


2.2. Temperature structure and heating terms 

The three-dimensional dust continuum radiative transfer code, 
RADMC-3D 1 is used to calculate the dust temperature structure. 

1 http://www.ita.uni-heidelberg.de/~dul1emond/ 
software/radmc-3d 


A v = 3 



Fig. 1 . Schematic showing the calculation of the thermal structure of the 
disk by combining both the RADMC3D Monte Carlo simulation and 
the diffusion equation. The Monte Carlo simulation calculates the tem¬ 
perature structure due to the irradiation (T in ) from a central star while 
the diffusion equation is used to solve the temperature structure at high 
optical depths (tr > 1 where tr is the Rosseland mean optical depth) 
as indicated by the red line. This high optical depth region starts typ¬ 
ically below the disk surface. The outflow wall, warm inner envelope 
(Tdust > 100 K and r < 500 AU) and photodissociation region are indi¬ 
cated. 



Fig. 2. Total luminosity (L tot = L* + L disk + L m ) as function of accretion 
rates for different stellar luminosities. 


A central star with temperature of 4000 K (a typical inferred ef¬ 
fective temperature of Class I protostars; White & Hillenbrand 
2004; Nisini et al. 2005) characterized by L* 1, 5 and 15 L Q is 
adopted. An accurate dust temperature structure of the disk is 
crucial in determining the location where the various volatiles 
thermally desorb from the grain. This is computationally chal¬ 
lenging for a massive optically thick disk being modelled here 
(see Min et al. 2009). Thus, we have separated the dust tempera¬ 
ture calculations due to the central star irradiation (passive) from 
the viscous heating treatment (see Fig. 1). The former is deter- 
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mined by RADMC3D considering a black body central star as 
listed in Table 1. 

An actively accreting disk provides additional heating (Ldi s k) 
from the loss of mechanical energy as the gas is viscously 
transported inward. The steady state accretion rate is typi¬ 
cally between 10~ 5 - 1(R 7 M 0 yr _1 (Hueso & Guillot 2005). 
However, episodic accretion events such as those simulated by 
Vorobyov (2009) can have transient spikes with an accretion 
rate up to 10 4 M 0 yr 1 . Thus, stellar accretion rates between 
10 8 —10 4 M 0 yr~' are adopted following the a disk formalism 
(Shakura & Sunyaev 1973; Williams & Cieza 2011). The vis¬ 
cous heating rate per volume is given by 

9 7 

Qv isc — ^Pdisk C’G 5 (4) 

where u = ac s II is the a-dependent turbulent viscosity parame¬ 
ter, c s the sound speed of the gas, and O the Keplerian angular 
velocity. At steady state, the viscosity is related to the disk mass 
and the accretion rate through (Lodato 2008): 

M = 3nvL. (5) 


To explore the degree of heating, the viscosity term v is varied 
for the explored accretion rates at a fixed disk mass (oc 2). The 
effective temperature (7\,j sc ) associated with the energy released 
at the inner radius assuming a hydrostatic disk is 


f' 


c sbL 4 (R) = Qvisc (R) dz - 


3 GM* 


-M 


1 - J — 


R * 

R 


(6) 


This is obtained by integrating the viscous heating terms verti¬ 
cally at all radii. Furthermore, 7\ isc is the effective temperature 
of the disk at the optically thin photosphere without the addition 
of stellar irradiation. However, the midplane temperature of an 
active disk is proportional to the vertical optical depth (Hubeny 
1990): 7" 4 ~ ^ R 2 gas r 4 isc with kr the Rosseland mean opacity 
which results in higher midplane temperatures. Such a method 
is similar to that of Kennedy & Kenyon (2008) and Hueso & 
Guillot (2005) in the optically thick regime. Note that the heat¬ 
ing from an accreting disk is caused by the dissipation of energy 
from both gas and dust. To account for the irradiation from the 
accreting disk, = f ncrT v i sc (R) 4 RdR is added to the central 
luminosity L* by determining its blackbody spectrum at T vlx at 
all radii of the disk. For the case of high accretion rates, it is ex¬ 
pected that the accretion proceeds onto the star since a gaseous 
disk is present within the dust sublimation radius. The total lu¬ 
minosity from the inner disk is L m = j R ncrT 4 isc RdR. Thus, 
the final irradiating central source L tot is the combined heating 
from the central star, the dusty disk, and the inner gaseous disk: 
Ltot = L* + Ldisk + L; n . Figure 2 shows the relation between L tot 
and M for different L* values. 

The following steps are taken to calculate the dust tempera¬ 
ture of an accreting embedded disk. 


- Monte Carlo dust continuum radiative transfer is used to 
simulate the photon propagation to determine the passively 
heated dust temperature structure (7’ ln ) due to total lumi¬ 
nosity (L tot ). The dust opacities of Crapsi et al. (2008) are 
adopted and are composed of a distribution of ice coated sil¬ 
icates and graphite grains. The removal of ices from the grain 
at Ldust > 100 K does not strongly alter the dust temperature 
structure. 


- Viscous heating is added to the region of the disk 
where the Rosseland mean optical depth r R > 1. 

This is done by fixing the midplane temperatures to 

Tmid = (^K R 2g as T 4 sc + L 4 r ) 1/4 . The vertical dust tempera¬ 
ture structure is calculated using the diffusion approxima¬ 
tion VDVL 4 = 0 bounded by T nr at the surface as ob¬ 
tained from the previous step and L m id at the midplane where 
D - Opdiisi'u-i ) • The diffusion is performed only within the 
r R > 1 regions. This calculation is repeated a few times such 
that it converges. The convergence is obtained when there is 
no change in the the r R > 1 region. 

The addition of the viscous heating can increase the dust temper¬ 
atures to > 2000 K while the typical dust vaporization tempera¬ 
ture is ~ 1500 K. Therefore, we have used an upper limit to the 
dust temperature of 1500 K. Gas opacities in the inner disk are 
not taken into account, which will affect the exact temperature in 
that region. This does not change the location of the snowlines 
since they are defined by dust temperatures 7’ ( | L1M < 160 K. 

The snowlines of protoplanetary disks without an envelope 
were obtained and compared with Min et al. (2011) to verify our 
approach (see Fig. A.l). Using this formulation, the differences 
in predicted water snowlines are typically within 2 AU at low 
accretion rates and less at high accretion rates. 

2.3. Molecular abundances 

The aim of this paper is to calculate the 2D snowlines or snow- 
surfaces for CO, CO 2 , and H 2 O in embedded disks. The region 
in which these volatiles freeze-out onto grains depends on the 
temperature and density structure (e.g., Meijerink et al. 2009). 
At steady state, the rate at which the molecule is adsorbed on 
the grain is balanced by the thermal desorption rate. We adopt 
dust number density Hd ust = 10~ 12 «h (Visser et al. 2009) with 
6: dust = 0.1 pm as the effective grain size. The number density of 
solids of species X (see also Fig. B.l) is simply given by 

nice = »dust7rfl 2 ust (3k B L gas /OT X ) 1/- ^ 

^gas D exp (—Lb/7dust) ^ 

where L gas gas temperature, 7d ust dust temperature, and mx is the 
mass of species X. The first-order pre-exponential factor, vi, is 
calculated from the binding energy. Lb (Hasegawa et al. 1992; 
Walsh et al. 2010) 


v 1 = 


2V ss Lb 

n-nix 


( 8 ) 


with the number of binding sites per surface area, N ss , is taken 
to be 8 x 10 14 cnT 2 following Visser et al. (2011). A dimen¬ 
sionless factor £(«i ce ) is used to switch between zeroth order to 
first order desorption when the ice thickness is less than a mono- 
layer. The properties for each molecule are given in Table 2 along 
with the calculated pre-exponential factor vi. These binding en¬ 
ergies assume pure ices. The typical timescales for adsorption 
are 4 - 6 x 10 3 year at temperature of 50 K and number densities 
(«h 2 ) °f 10 6 cm 4 . This implies that the steady state assumption 
is not valid at lower densities present in the large-scale enve¬ 
lope (r > 1000 AU) where the freeze-out timescales becomes 
longer than the lifetime of the core (e.g., Jprgensen et al. 2005). 
On the other hand, this paper focuses on the gas phase abun¬ 
dances at small-scales r < 100 AU where the number densities 
are «h 2 > 10 6 cm -3 . Photodesorption can be ignored at such 
high densities and within the disk, but it may be important at the 
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Table 2. Molecular parameters to calculate R its as tabulated in Burke 
& Brown (2010). 


Molecule 

n a 

[s’ 1 ] 

Eb 

[K] 

Refs. 

CO 

6.4 x 10 11 

855 

Bisschop et al. 2006 

C0 2 

8.2 x 10 11 

2400 

Galvez et al. 2007 

h 2 o 

2.1 x 10 12 

5773 

Fraser et al. 2001 


Notes. (al First order desorption pre-exponential factor calculated from 
the binding energies. 

disk’s surface and along the outflow cavity wall. Since the water 
vapor will also be rapidly photodissociated in those locations, 
these regions are not major water reservoirs (see Fig. 1). To ap¬ 
proximate this region, water is assumed not to be present within 
a region that is characterized by Ay < 3 or V H < 6x 10 21 cm~ 2 . A 
more detailed gas phase abundance structure through a chemical 
network will be explored in the future. 

3. Results 

3.1. Thermal structure of an actively accreting embedded 
disk 




Fig. 3. Midplane radial ( bottom) and vertical temperatures (top) at 1, 
5, and 25 AU for an embedded 0.1 M 0 disk and an accretion rate 
of 10~ 5 M 0 yr~' irradiated by L* = 1 L e . The red dashed ( bottom ) 
and solid (top) lines indicate the thermal structure of a passively ir¬ 
radiated disk (T m ) calculated by the Monte Carlo simulation. The 
blue solid lines indicate the temperatures including the viscous heating 

(Tefr = (^KR£ gas 7A sc + T^j ). The viscous temperature as calculated 
from Eq. 6 is indicated by the dotted blue lines. 

The locations where various molecules can thermally desorb 
from the dust grains depend on the temperature structure of the 
disk. Irradiated disks have a warm upper layer with a cooler mid¬ 
plane. Figure 3 (red dashed lines) shows the midplane ( bottom ) 
and vertical temperature at a number of radii (top) for an em¬ 
bedded disk passively irradiated by a 1 L Q central source. Here, 



z/r 


Fig. 5. Comparison of the vertical temperature structure at 1, 5, 25, and 
150 AU between a protoplanetary disk (red) and embedded disks (blue 
and black). The difference between the two embedded disk models is 
the centrifugal radius R ccn = 50 AU (black) versus 200 AU (blue). The 
envelope mass is 1 M e for the embedded disk model. 


we present the results for a 0.1 M 0 disk embedded in a 1 M 0 
envelope whose R cen is 200 AU (Fig. 3 top). These canonical 
parameters are highlighted in Table 1. 

The dust temperatures of an embedded actively accreting 
disk are indicated by the blue lines in Fig. 3. The dotted line 
shows the viscous temperature (7' v j sc ) as expected at the photo¬ 
sphere (optically thin) while the solid blue lines show the effec¬ 
tive dust temperature corrected for the optical depth and passive 
irradiation. As previously found, the addition of viscous heating 
can raise the temperatures in the inner few AU to > 1000 K (e.g., 
Calvet et al. 1991; D’Alessio et al. 1997; Davis 2005). The disk 
temperature is above the water sublimation temperature out to 
R ~ 30 AU for an accretion rate of 10 5 M 0 yr 1 . Due to the 
R 3 dependence of r v ; sc (see Eq. 6), the viscous heating is dom¬ 
inant in the inner few AU as indicated in Fig. 3 (bottom). Conse¬ 
quently, the passive irradiation due to the central luminosity (L*) 
dominates the temperatures along the disk’s photosphere and the 
outflow cavity wall while the viscous dissipation dominates the 
heating deep within the disk. These effects can be seen in the 2D 
dust temperature structure in the inner ~ 40 AU shown in Fig. 4 
for three different accretion rates. For M > 10 5 M 0 yr _1 , panels 
show the vertical temperature inversions. 

The dust temperature structure is also compared with a disk 
with and without an envelope whose R cen is 50 AU in Fig. 5. 
The difference is small in the inner disk, and lies primarily at 
large radii (R >1 AU) where the midplane temperature structure 
(7 m id ce £) is weakly affected by the adopted envelope model. 
The temperature structure of the embedded disk model depends 
slightly on the adopted envelope model. Most importantly, at 
R — 5 and 150 AU, the embedded disk model with R cen = 50 
AU is warmer than the other two models. This is mainly due to 
the mass distribution in the inner 100 AU. Overall, however, the 
differences are small (see also D’Alessio et al. 1997) and, for 
our purposes, it is sufficient to fix the envelope mass to 1 M 0 
to assess the overall temperature stmcture of an embedded ac¬ 
tively accreting disk. The effects are further discussed in §3.4. 
However, the difference at larger radii is sufficient to affect the 
dominant phase of CO 2 and CO. 
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Fig. 4. Dust temperature structure in the inner 44 AU for three different accretion rates for a 0.1 M a disk embedded in a 1 M a envelope. A 1 
L q central heating source is adopted for these models. The disk is oriented horizontally while the outflow cavity is oriented vertically. The three 
different lines indicate 160 (dashed-dot), 100 (dashed), and 50 (solid white) K contours, which are important for the water snowlines (see Fig. B. 1). 


3.2. Water snowline 



Fig. 6. Midplane water snowlines as function of accretion rates and disk 
mass for a 200 AU embedded disk with f? ccn = 200 AU. The different 
lines indicate the snowlines dependence on disk mass at a fixed stellar 
luminosity and envelope mass. 

Using the obtained dust temperatures, the gas and ice number 
densities are calculated at each cell. To determine the snowline, 
the total available water mass was computed by adopting a water 
abundance of 10 4 with respect to Hi. The total available mass is 
then multiplied by the gas fraction to determine the total water- 
vapor mass. A minimum gas water abundance of 10 with re¬ 
spect to FD has been used. The snowline is defined as the radius 
at which 50% of the total available water has frozen onto the 
grains (M gas /Mi ce = 0.5) (e.g., Min et al. 2011). This typically 
occurs at ~160 K in the high density regions (n H ~ 10 14 cm~ 3 , 
see dashed line in Fig. 4). 

Figure 6 presents the midplane snowline radius as a function 
of accretion rate. In the absence of accretion heating, the water 
snowline is located at ~ 3 AU. This does not strongly depend on 


the accretion rate until a value of M > 10 7 M 0 yr 1 is reached. 
The maximum water snowline is located at 55 AU for an accre¬ 
tion rate of 10 4 M 0 yr 1 . Although the effective disk midplane 
temperature depends on the disk mass, it does not strongly affect 
the water snowline location. It is located at only slightly smaller 
radius for a less massive disk (0.05 M 0 ) as indicated in Fig. 6. 

The steep increase of the water snowline at high accretion 
rates can be understood by comparing the stellar luminosity and 
the accretion luminosity. For the canonical values that are in¬ 
dicated in Table 1, the irradiating central luminosity is 1 L 0 . 
The accretion luminosity is estimated by integrating Eq. 6 ra¬ 
dially over the active disk between 0.1 and 200 AU, in this 
case, and is approximately L acc ~ 0.5 x GM stm M/R[ n . Thus, 
the accretion luminosity is equal to that of the central star for 
M ~ 3 X 10 6 M q yr 1 . The addition of the accretion onto the 
star due to the gaseous inner disk provides an additional ~ 10 L 0 
luminosity at the given accretion rate. Thus, the accretion lumi¬ 
nosity starts to contribute to the heating at M > 10~ 7 M 0 yr~' 
for the adopted parameters (see Fig. 2). 

The observable water emission depends on the water vapor 
column density. The water vapor column extends further than 
the snowline due to the vertical gradient in water vapor abun¬ 
dance. The available water is rapidly frozen out onto dust grains 
beyond the snowline. To determine whether the location of the 
snowline is within the disk or not, the hydrostatic disk surface 
is determined through H = c s /Qk where c s = pniu is 

the sound speed and Ok the Keplerian angular frequency. This 
is the approximated regime where the gas should be in Keple¬ 
rian motion. For a 1 L 0 central star, most of the water emission 
arises from the warm inner envelope (sometimes also called the 
‘hot core’) above the hydrostatic disk surface as shown by the 
yellow line in Fig. 7 if a constant water abundance of 10 4 with 
respect to Hi is adopted. The blue shaded region in Fig. 7 indi¬ 
cates where water vapor is dominant, some water vapor (< 50% 
in mass) is still present within the inner few AU of the red re¬ 
gions (< 30 AU). Water vapor is also present along the cavity 
walls as indicated by the two purple lines in Fig. 7, however 
these regions have low extinctions (A v < 3) which allows for the 
photodissociation of water. 
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Fig. 8. Left : Midplane CO 2 snowlines as a function of accretion rates and disk mass for a 200 AU disk. Right: Midplane CO snowlines as a function 
of accretion rates for the same disk. The different colors indicate the snowlines dependence on disk mass at a fixed stellar luminosity (1 L 0 ) and 
envelope mass (1 M 0 ). 


M = 10“ 5 M 0 yr —1 



R [AU] 


Fig. 7. Locations of gas phase volatiles in the 0.1 M a disk embedded 
in a 1 M 0 envelope. The different colors indicate the different regimes 
where various volatiles are found in the gas phase: blue (H 2 0, C0 2 , 
and CO), red (only C0 2 and CO), and green (only CO). The accretion 
rate is indicated at the top. The arrows are indicating three different 
regions as indicated in Fig. 1: water photo-dissociation region defined 
by A v = 3, warm inner envelope, and the disk surface (yellow line 
assuming T gas = T dust ). 


3.3. CO and C0 2 snowlines 

Pure CO 2 and CO ices thermally desorb over a narrow range of 
dust temperatures between 40-80 K and 15-30 K, respectively, 
depending on the density (see Fig. B.l in the appendix). Due to 
their lower binding energies relative to water, CO and CO 2 are 
in the gas phase within a large part of the embedded disk. The 
CO 2 snowline is between 20-250 AU for 10 9 —10 4 M 0 yr 1 


accretion rates (see Fig. 8). Thus, the entire disk including the 
midplane lacks CO 2 ice for highly acccreting embedded disks 
(M ~ 10 -4 M g yr -1 ). The steep rise of the CO 2 snowline at high 
accretion rates is similar to that of water as shown in Fig. 6. 

Since the adopted binding energy of CO to the dust grain is 
the least with 855 K, CO largely remains in the gas phase within 
the disk for L* > 1 L Q independent of the accretion rate (see 
Fig. 8 right). At low accretion rates, the snowline is located at 
~ 150 AU at the midplane indicating the presence of CO ice be¬ 
tween 150 to 200 AU within the disk. However, the bulk of CO 
within the disk remains in the gas phase as shown in Fig. 7. A 
smaller disk during the embedded phase would lead to stronger 
envelope irradiation (e.g., D’Alessio et al. 1997), and, conse¬ 
quently, CO is not frozen out within the disk. A sufficiently large 
and massive disk (M^k > 0.2 M G irradiated by a 1 L Q star) 
could contain a larger fraction of CO ice at large radii. This is 
simply due to the increase of optical depth and, thus, lower dust 
temperatures at large radii. 

3.4. Dependence on stellar ; disk , and envelope properties 

Previous sections presented the vapor content for a 0.1 M 0 em¬ 
bedded disk surrounded by a 1 M 0 envelope with R cen = 200 AU 
that is being irradiated by a 1 L 0 central star (7’* = 4000 K, 
R+ =2.1 R q ) as highlighted in Table 1. As noted in Section 2.2, 
the effective midplane temperature (T^\) depends on the disk 
mass but an increasing disk mass leads to only a small change 
to the water snowline as shown in Fig. 6. Figure 7 shows that 
the water vapor can be either within the disk or within the warm 
inner envelope. This section explores how the vapor content de¬ 
pends on disk radius and central luminosity for fixed envelope 
and disk (0.1 M 0 ) masses. 

Figure 9 presents the locations of water snowline as a func¬ 
tion of accretion rate, luminosity, and disk radius. For a fixed 
disk structure, a factor of 5 increase in L* shifts the midplane 
water snowline by a factor of ~2 from ~4 AU (1 L Q ) to ~10 AU 
( 10 L e ) at low accretion rates (M < 10 -7 M 0 yr -1 ). A factor of 
15 increase in L* yields a factor of ~ 4 increase in the location 
at the water snowline. The stellar luminosity L* is defined by the 
stellar radius while the accretion dependent total luminosity L tot 
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Fig. 9. Midplane water snowline as a function of stellar luminosity, accretion rate, and disk radius. The envelope mass is fixed at 1 M e with a disk 
mass of 0.1 M e and f? ccn = 200 AU. The different panels show the water snowline as a function of stellar luminosity: 1 L 0 (left), 5 L 0 (center), and 
15 L Q (right). The different colors indicate the disk radius: 200 AU (blue), 100 AU (green), and 50 AU (red). 


that is irradiating the system is shown in Fig. 2. For high accre¬ 
tion rates (M > 10 s M 0 yr -1 ), the water snowline decreases as 
L* increases. This can be explained due to the fact that L tot de¬ 
creases as L* increases at high accretions because /?* increases, 
which in turn decreases the luminosity that is provided by the 
inner gaseous disk. Thus, if the stellar accretion rates are low 
(L* > Lace), the stellar luminosity is the parameter to be var¬ 
ied in order to increase the water snowline. However, the water 
snowline can be up to 50-60 AU if the stellar luminosity is low 
(small stellar radius) and the stellar accretion rate is high. 

The CO 2 midplane snowline shows similar behavior as the 
water snowline (Fig. B.2) as function of stellar luminosity and 
accretion rates. It is in the gas phase at R > 40 AU for L* > 5 L Q , 
which is a factor of 2 increase in terms of snowline location with 
respect to the canonical value of L* = 1 L G . As CO is already 
largely in the gas phase for the canonical values, a small increase 
of the stellar luminosity results in the lack of a CO snowline 
within the disk (Fig. B.3). 

Another parameter that affects the total energy input into the 
embedded system is the disk radius. The total accretion lumi¬ 
nosity that is being added depends on the disk radius, which is 
the region where the disk’s accretion energy is being integrated. 
A decreasing disk’s radius leads to a decreasing accretion lu¬ 
minosity. This only affects the midplane water snowline in the 
case of L* = 1 L 0 ( see Fig. 9 ). In this case, the difference is 
seen only for high stellar accretion rates where the accretion lu¬ 
minosity (Ldi s k + L; n ) starts to be the dominant heating source. 
The midplane water snowline decreases from ~ 55 AU to ~ 38 
AU as the disk radius decreases from 200 AU to 50 AU. Thus, an 
increasing stellar luminosity (L*) has a greater effect on the mid¬ 
plane water snowline. This leads to an increase of water vapor in 
the warm inner envelope (see Fig. 7) for low stellar accretion 
rates (M < 10 -5 M 0 yr -1 ). However, if L* > 1 L 0 , the disk ra¬ 
dius does not affect the location of the midplane water snowline 
even for high accretion rates (M > 10 -5 M Q yr -1 ) since the to¬ 
tal luminosity is dominated by the stellar luminosity. The CO 2 
snowline, on the other hand, shows a stronger dependence on 
the disk radius than that of water. For the high stellar luminos¬ 
ity case L* = 15 L 0 , the CO 2 snowline increases from 40 AU 
to 80 AU as the disk radius (/?disk) increases from 50 AU to 200 


AU. This can be explained by the density distribution in the inner 
region. For a small disk, the dust temperature is low enough at 
large radii such that CO 2 can exist in the solid phase. However, 
as the size of the disk increases, the overall density at similar 
radii decreases which results in lower optical depth and, as con¬ 
sequence, higher dust temperatures at large radii. The disk radius 
parameter seems to be more important for the CO 2 snowline than 
for the water snowline. 

The main parameter for the envelope in our setup is the cen¬ 
trifugal radius. In the case of the canonical value of R cen = 200 
AU, the mass of the disk can be distributed between 50, 100, 
or 200 AU. The disk’s radius /Lusk is where the disk ends and 
it can be smaller than R cen . Here, the effects of a flattened en¬ 
velope as defined by the centrifugal radius R cen are presented. 
A smaller centrifugal radius results in a hotter inner disk (see 
Fig. 5). However, the affected regions correspond to the regions 
within the disk where dust temperatures are already > 100 K. 
Thus, we find that there is no change in the water snowline for 
a more flattened envelope stmeture (R cea = 50 AU instead of 
200 AU). More volatile species such as CO 2 and CO are affected 
by the flattened envelope structure, however. The CO 2 snowline 
is confined within the inner 40 AU for M < 10 -5 M 0 yr -1 . At 
higher accretion rates, the combination of a warm disk and low 
envelope densities at R > 50 AU for the R cen = 50 AU model 
results in a warm envelope where CO 2 is largely in the gas phase 
up to 400 - 500 AU along the midplane while it is ~ 200 AU in 
the canonical model. In the case of CO, its snowline is already 
at 400 AU for an R cen = 50 AU envelope model at low accre¬ 
tion rates instead of ~ 200-300 AU in the case of R cen = 200 
AU envelope. Thus, we find that the inner envelope and disk 
physical stmeture to be important in determining the location of 
snowlines of various volatiles. Yet, these parameters are largely 
unknown for deeply embedded young stellar objects. 

4. Discussion 

4. 1. Comparison with accreting protoplanetary disk models 

In the absence of an envelope, similar stellar parameters and an 
accretion rate of 10 -8 M 0 yr -1 lead to the midplane water snow¬ 
line located at ~ 1 AU (e.g., Sasselov & Lecar 2000). Lecar et al. 
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(2006) suggest that the snowline can move out by increasing the 
accretion rate, disk mass and the dust opacities. The differences 
between the derived snowlines in the literature are unlikely to be 
due to the differences in radiative transfer treatment which has 
been compared in Min et al. (2011). As an example, Garaud & 
Lin (2007) derived similar water snowlines to that of Min et al. 
(2011) (1 AU vs 2 AU) at similar accretion rates using different 
methods in deriving the dust temperatures and different opaci¬ 
ties. The latter model takes into account the detailed 2D vertical 
structure of the disk and dust sublimation front. As Min et al. 
(2011) show, most of the differences are due to the adopted dust 
opacities. Our value is ~ 1 AU larger than the values tabulated 
in Min et al. (2011) for similar parameters. Furthermore, most 
of the previous studies adopt the minimum mass solar nebula 
(MMSN) model where 2 oc R 1 5 instead of the R 1 2D para¬ 
metric disk structure used in this paper (see e.g., Andrews et al. 
2009). Previous snowlines or snow surfaces studies do not take 
the flattened envelope into account. 

Previous studies have found that the water snowlines in the 
more evolved disks decrease with decreasing accretion rates 
(e.g., Davis 2005; Garaud & Lin 2007; Min et al. 2011). This 
has been reproduced in Fig. A.l. However, if the disk structure 
is related to its stellar accretion rate such as adopted in Garaud 
& Lin (2007), the water snowline increases to a constant value 
at ~ 2 AU for very low accretion rates (M < 10 -10 M 0 yr -1 ). 
This is due to the disk becoming more optically thin, which in 
turn allows for stellar photons to penetrate deeper radially into 
the disk. In this paper, the disk structure is fixed while the stellar 
accretion rate is varied. However, the envelope provides a blan¬ 
ket where the disk can stay relatively warm. Thus, even though 
the disk structure is fixed, the snowline stays at a fixed radius for 
decreasing stellar accretion rates. 

Taking the envelope into account, for the same accretion rate 
of 10 8 M q yr 1 , our water snowline is at ~ 3 AU for embedded 
disks compared to 1-2 AU without an envelope. The presence of 
the envelope does not strongly affect the midplane temperature at 
regions close 100 K as shown in Fig. 5. Thus, it does not modify 
the water snowline significantly. However it does affect the CO 2 
and CO snowlines since the 40 K region shifts inward under the 
presence of an envelope depending on the centrifugal radius. 

4.2. Caveats 

The effect of convection in the vertical direction is not included 
in this study. It is typically found to be important in the case 
of high accretion rates (M > 10 -6 M 0 yr -1 , D’Alessio et al. 
1998, Min et al. 2011). Convection is found to cool the midplane 
temperatures at 7'<| USl > 500 K. Thus, this should not affect the 
water (100-160 K), C0 2 (~ 50 K) and CO (~ 20 K) snowlines 
as their sublimation temperatures are much lower than 500 K. 

Under the assumption of the steady state accretion disk 
model, the typical values of a are found to be > 1 for the 
case of an accretion rate of M — 10 4 M 0 yr -1 . This is signif¬ 
icantly larger than that expected from magnetorotational insta¬ 
bility (MRI) driven accretion in a MHD disk (a = 0.01, Balbus 
& Hawley 1998). King et al. (2007) indicated that a could be 
~ 0.4 for thin disks around compact systems. Such large a val¬ 
ues in our models are obtained since we have fixed the density 
distribution at a high accretion rate (M > 10 5 M 0 yr -1 ). In real¬ 
ity, such a disk is unphysical. However, this is only encountered 
for the highest accretion rate of 10 4 M 0 yr -1 that occurs for 
a very short time. In principle, this can be avoided by iterating 
the disk structure to ensure that the hydrostatic conditions are 
satisfied. However, this paper explores the effect of the adopted 


accretion rates and disk masses for a given fixed set of parame¬ 
ters. Self-consistent models including gas opacities and multiple 
dust species should be explored in the future. 


4.3. Comparison with observations 

The extent of the observed spatially resolved water emission 
toward low-mass embedded YSOs is between 25-90 AU (J0r- 
gensen & van Dishoeck 2010; Persson et al. 2012, 2014). The 
water snowline is confined to within the inner 40 AU for the 
cases of the envelope models whose centrifugal radius is 50 AU. 
In such models, we have chosen that the disk does not extend 
beyond 50 AU. The water snowline is only increased in cases 
where the centrifugal radius of the envelope is larger than 50 
AU, in which case the water snowline can extend up to 55 AU 
if the stellar luminosity is low (L* = 1 L e ). Thus, our results 
indicate that most of the water emission at R < 50 AU would 
originate from the combined disk and warm inner envelope (see 
Fig. 7). 

Through the combined modelling of the continuum spec¬ 
tral energy distribution and submm continuum images, the enve¬ 
lope masses toward the observed embedded sources (NGC1333 
IRAS2A, IRAS4A, and IRAS4B) are found to be between 3-5 
M 0 (Kristensen et al. 2012), somewhat larger than our canoni¬ 
cal value. On the other hand, all three sources are binaries (e.g., 
Jprgensen et al. 2004, 2007) so it is likely that the effective en¬ 
velope mass that affects the disk physical structure is lower than 
the above values. Therefore, the envelope masses have been kept 
fixed at 1 M 0 . 

The amount of water vapor in the embedded disk models 
can be compared to the observed optically thin millimeter wa¬ 
ter emission toward deeply embedded YSOs. Thermalized H 18 0 
emission at 203.4 GHz (3 i, 3 - 2 2 ,o) is calculated using the fol¬ 
lowing equation 


Iy ~ B v (T ex )( 1-0, 

_ _ N A u ic 3 hv/kRTn 

Q(T ex )AV 8ttv 3 


(9) 

( 10 ) 


where A u \ is the Einstein A coefficient of the 203 GHz transition 
is 4.5 x 10 -7 s _1 , r ex is taken to be 150 K, E u = 204 K is the 
upper energy level, and Q is the temperature dependent partition 
function adopted form the HITRAN database (Rothman et al. 
2009). The beam averaged column density (N) is taken from 
within 50 AU radius with l6 0: 18 0 isotopic ratio of 540 (Wil¬ 
son & Rood 1994). The line profile is assumed to be Gaussian 
with a FWHM (A u) = 1^1 km s 1 as observed toward the three 
embedded sources reported by Persson et al. (2014). 

Figure 10 (left) presents the mean H ls O gas column den¬ 
sity of the simulated embedded disks as a function of accre¬ 
tion rate for different disk masses and stellar luminosities. The 
mean column density rises sharply for M > 10 7 M 0 yr -1 as ex¬ 
pected for a 1 L Q stellar luminosity while the sharp rise occurs 
at M > 10 7 M q for a 5 L 0 irradiating source. The typical line 
optical depth at the line center (u = 0 km s -1 ) can be t v= 0 > 1. 
The typical average column densities of water between 50 AU 
to 1000 AU are > 4 orders of magnitude lower than the values 
within 50 AU as presented in Fig. 10. Thus, the water present 
within 50 AU dominates the total water mass for the adopted 
models. 

The expected integrated H 18 0 3 1 . 3 —2 2 ,o (203.4 GHz) line flux 
densities for the embedded disk models are also shown in Fig. 10 
(right). The integrated line flux densities (J S Y dv) are calculated 
with a Gaussian line profile with Av — 1 km s -1 (1.06 S v= 0 x Av 
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Fig. 10. Left: Beam averaged water (H 18 0) column densities within 50 AU radius as a function of accretion rates. The different colors indicate 
the disk mass dependence for R disk = 200 AU models: M d = 0.05 M a (red triangles), 0.2 (blue circles), and 0.5 (green squares). The different 
lines show the luminosity dependence: L* = 1 L 0 (solid) and 5 L 0 (dashed). Center: Integrated line flux densities arising from the disk as defined 
by the yellow line in Fig. 7. The observed integrated H^ 8 0 line is indicated by the gray lines for the three different embedded YSOs in Persson 
et al. (2014). Right: Integrated line flux densities accounting for the entire water vapor mass assuming a Gaussian linewidth of 1 km s -1 which is 
appropriate for NGC1333-IRAS4B. 


where S v= o is the peak flux density at line center) as observed 
toward NGC1333-IRAS4B (Persson et al. 2014). The observed 
line widths toward IRAS2A and IRAS4A are 4 and 3 km s _I , re¬ 
spectively. Increasing A u slightly increases the model integrated 
line flux densities because the line is slightly optically thick. 

The predicted integrated line flux densities from actively 
accreting embedded disk models are consistent with that ob¬ 
served toward NGC1333 IRAS4B if Ma < 0.2 M 0 with M > 
6 x 10- 7 M 0 yr 1 for L* = 1 L 0 . Jprgensen et al. (2009) esti¬ 
mated from the continuum that the disk mass is ~ 0.2 M G , which 
is consistent with our results. Decreasing the disk mass implies 
a higher accretion rate is required to reproduce the observed wa¬ 
ter emission. A higher accretion rate (~ 1 x 10 7 M 0 yr 1 ) with 
lower disk mass is also consistent with the observed extent of 
water emission at ~ 20 AU. A similar conclusion is reached for 
the case of IRAS4A whose disk is estimated to be ~ 0.5 M 0 
(Jprgensen et al. 2009). Since the bolometric luminosities toward 
NGC1333 IRAS4A and IRAS4B are a factor of 4-9 higher than 
the adopted central luminosity (L*), a model with a factor of 2 
lower mass disk can also reproduce the observed flux density 
with accretion rates ~ 10 s M Q yr -1 . The total luminosity gener¬ 
ated by such embedded systems would be similar to the observed 
bolometric luminosities. With such a lower mass disk, the warm 
inner envelope is the dominant source of emission instead of the 
Keplerian disk. This implies that it would be difficult to observe 
the Keplerian motion of the disk from H ls O lines. 

IRAS2A shows stronger integrated line emission (~ 1 Jy km 
s 1 ) relative to the other 2 sources and extending up to 90 AU 
radius. As pointed out above, a disk cannot be responsible for 
any warm water emission at R > 50 AU. Furthermore, Brinch 
et al. (2009) estimated that the disk is at most ~ 0.05 M 0 , which 
lowers the expected disk’s contribution to the observed water 
emission (see also Jprgerisen et al. 2009 and Maret et al. 2014). 
There are two possible scenarios that can explain such a high wa¬ 
ter emission: a more luminous star (L* > 10 L 0 ) or a high stellar 
accretion rate (~ 10 s M 0 yr -1 ). Both solutions result in total lu¬ 
minosity L u „ comparable to the observed bolometric luminosity 


of IRAS2A (L ho i = 35.7 L 0 , Karska et al. 2013). In either case, 
most of the water emission arises from the warm inner envelope 
with negligible disk contribution. 

The best cases for an actively accreting embedded disk sce¬ 
nario are IRAS4A and IRAS4B. More detailed modelling of 
their physical structure at < 100 AU radius is required to fur¬ 
ther constrain the relation between infall (envelope to disk) and 
stellar accretion rate (disk to star). Mottramet al. (2013) suggests 
that the infall rate from the large-scale envelope toward IRAS4A 
at 1000 AU to be ~ 10 -4 M 0 yr -1 . For the cases of disk masses 
between 0.2-0.5 M 0 , the disk must also process the material at 
similar rate (within a factor of 10), which is consistent with the 
models presented here. 

Observational perspectives Future spatially and spectrally re¬ 
solved ALMA observations of optically thin water (H 18 0 ) and 
deuterated water (HDO) lines may be able to differentiate these 
models. The spatially resolved data can provide the extent of the 
water emission. The comparison between the source’s bolomet¬ 
ric luminosity and the extent of the water emission can provide 
limits on the stellar accretion rate. A low bolometric luminosity 
source cannot have a water emission originating from the disk 
further than a few AU. Furthermore, the velocity information 
within the inner 50 AU radius can help in differentiating between 
the disk and the infalling envelope (see e.g., Harsono et al. 2015). 

4.4. Connecting to the young solar nebula 

It is thought that the disk out of which the solar system formed 
(the solar nebula) was initially hot enough to vaporize all ma¬ 
terial inherited from the collapsing cloud into atoms, then turn 
them into solids according to a condensation sequence as the 
nebula cools (Lewis 1974; Grossman & Larimer 1974). Evi¬ 
dences for energetic processing and subsequent condensation 
comes from meteoritic data collected in the inner solar system 
(see Kerridge & Matthews 1988; Scott 2007; Apai & Lauretta 
2010, for reviews). After the refractory phases have formed. 
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more volatile species such as ices can also condense at larger dis¬ 
tances from the young Sun, with ice composition depending on 
the temperature and pressure as well as the elemental abundance 
ratios (e.g., Lunine et al. 1991; Owen & Bar-Nun 1993; Mousis 
et al. 2009; Pontoppidan et al. 2014). Using cooling curves ap¬ 
propriate for the solar nebula, models show that various ices in¬ 
cluding PUO, CH 3 OH, N'Hi and CO 2 ice form by condensation 
out of warm gas out to at least 20 AU (Mousis et al. 2012; Mar- 
boeuf et al. 2014). This condensation process is important to ex¬ 
plain the composition of volatiles in the atmospheres of solar- 
system giant planets and comets. 

Where and when does this heating and condensation in the 
disk actually occur? The above processes are usually discussed 
in the context of protoplanetary disks where the envelope has 
dissipated and where the water snowline is eventually at a few 
AU once the disk has cooled. However, neither observations nor 
models of disks around solar-mass T-Tauri stars show any evi¬ 
dence that disks are as warm as required above to have conden¬ 
sation happening out to 20 AU (Beckwith et al. 1990; D’Alessio 
et al. 1998; Dullemond et al. 2007). Furthermore, observations 
suggest that significant grain growth has occurred (Testi et al. 
2014) and disk inhomogeneities due to protoplanets are present 
in such evolved disks (e.g., Quanz et al. 2013). This further vi¬ 
olates the assumption that the physical structure of an evolved 
disk can be used to study the young solar nebula. The results 
from this paper indicate that instead such a hot disk may be 
more common in the early stages of disk formation during the 
embedded phase (M env > The disk is then expected 

to have high accretion rates providing the necessary additional 
heating. Furthermore, the high accretion rates necessary to push 
the midplane snowlines to 20 AU and beyond tend to occur 
for only a short time (e.g., Vorobyov 2009). Most volatiles in¬ 
cluding H 2 O, CH 3 OH, CO 2 and CO presented here will then be 
in the gas-phase within this 20-30 AU radius through sublima¬ 
tion of ices. Shortly after, the volatiles re-condense as the accre¬ 
tion rate decreases, with the process controlled by the freeze- 
out timescale. The inclusion of the gas motions (radially and 
vertically) and kinetics does not allow for instantaneous re¬ 
condensation, thus some gas-phase volatiles should remain near 
the snowline (Lewis & Prinn 1980; Ciesla & Cuzzi 2006). 

A related question is whether the volatiles that are incorpo¬ 
rated into planetesimals and eventually planets and icy bodies 
in the critical 5-30 AU zone are then inherited or reset during 
the disk formation process. The material is defined to be inher¬ 
ited if the simple and complex ices that have been built up in 
the envelope survive the voyage to the planetary and cometary 
forming zones. Evolutionary models by Visser et al. (2009) and 
Visser et al. (2011) suggest that strongly bound ices such as H 2 O 
are largely pristine interstellar ices in the outer disk, whereas 
more volatile species can have sublimated, recondensed and re¬ 
processed several times on their way to the inner disk. However, 
these models have not included internal viscous dissipation as 
an additional heating source. A warmer early disk can facilitate 
ice evolution and formation of complex organic molecules in the 
temperature regime between 20-40 K (Garrod et al. 2008; No¬ 
mura et al. 2009; Walsh et al. 2014; Drozdovskaya et al. 2014). 
Our results including the accretion heating indicate that signifi¬ 
cant ice re-processing may occur out to larger radii than thought 
before, well in the comet formation zone. Such a reset scenario 
is supported by the meteoritic data from the inner few AU, but 
the fraction of material that is reset at larger radii is still un¬ 
known (e.g. Pontoppidan et al. 2014). If some of the icy grains 
are incorporated into km-sized bodies (planetesimals) immedi¬ 


ately following disk formation, a larger fraction of the material 
could still remain pristine. 

5. Summary and conclusions 

Two-dimensional embedded disk models have been presented to 
investigate the location of the snowlines of H 2 O, CO 2 and CO. 
The dust temperature structure is calculated using the 3D dust 
radiative transfer code RADMC3D with a central heating source 
(L* +L v isc + Ein)- An additional heating term from an actively ac¬ 
creting disk has been added through a diffusion approximation. 
The main parameters that are explored in this paper are L*, cen¬ 
trifugal radius, disk radius, disk mass, and stellar accretion rate. 
In addition, the extent of the snowline and optically thin water 
emission are compared with observations toward three deeply 
embedded low-mass YSOs. The following lists the main results 
of this paper. 

- The presence of an envelope serves as a blanket for the 
disk such that the dust temperature within the inner few AU 
stays relatively warm for the adopted disk parameters. The 
adopted centrifugal radius affects the temperature structure 
of the disk since the effective midplane temperature depends 
on the optical depth. An envelope with a smaller centrifu¬ 
gal radius distributes more mass in the inner 100 AU than an 
envelope with a higher centrifugal radius. Such high concen¬ 
tration of mass in the inner 100 AU increases the effective 
Rosseland mean optical depth that an observer at the disk’s 
midplane sees vertically, which in turn increases the effective 
midplane temperature. While it does not drastically affect the 
?dust > 100 K regime, the presence of an envelope affects the 
location of 7 duM < 40 K, which is important for the CO and 
CO 2 snowlines. 

- The midplane water snowline in the models can extend up to 
a maximum of ~ 55 AU for a disk accreting at 10 4 M 0 yr 1 . 
The CO 2 snowline is located at >100 AU for the same accre¬ 
tion rate. Both H 2 O and CO 2 can remain in the solid phase at 
large radii (R > 100 AU) in the midplane within the bound¬ 
aries of an embedded hydrostatic disk. 

- CO is largely found to be in the gas phase within the em¬ 
bedded disk independent of accretion rate and disk prop¬ 
erties. Some CO could be frozen out at large radii in the 
midplane but only for a relatively massive embedded disk 
(Mi > 0.1 M 0 ). 

- The CO 2 and CO snowlines are affected by the exact struc¬ 
ture of the flattened envelope. A smaller centrifugal radius 
corresponding to a highly flattened envelope model increases 
the CO and CO 2 snow surfaces to larger radii. 

- The observed H' s O emission toward NGC1333 IRAS4A and 
IRAS4B is consistent with models for M > 10 6 M e yr -1 for 
the adopted disk masses of 0.5 M 0 and 0.2 M 0 , respectively. 
The observed size of water emission can also be explained 
with a model with a slightly higher accretion rate but with a 
lower disk mass (R < 50 AU). Toward IRAS2A, the models 
suggest a smaller disk contribution than toward IRAS4A and 
IRAS4B. Most of the optically thin water emission is emitted 
from the warm inner envelope (hot core) for that source ( see 
Fig. 7). 

- Significant chemical processing is expected to occur in the 
inner 100 AU region during the disk formation process in 
Stage 0. Midplane temperatures between 2(M10 K are ex¬ 
pected to be prevalent up to 100 AU radius, and 100 K out to 
30 AU, but only for the highest accretion rates found in the 
early phases of star formation. During these warm phases. 
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the chemistry inherited from the collapsing cloud could be 
reset and more complex molecules could form, unless the 
ices are sequestered early into planetesimals. 

In connection to early solar system formation, an important 
conclusion from our models is that a hot young solar nebula with 
water vapor out to 30 AU can only occur during the deeply em¬ 
bedded phase, not the T Tauri phase of our solar system. Our 
models also show that most of the observed optically thin wa¬ 
ter emission arises within the inner 50 AU radius. Future Ata¬ 
cama Large Millimeter/submillimeter Array (ALMA) observa¬ 
tions spatially resolving the inner 50 AU can test our embedded 
accreting disk models when coupled with better physical models 
for individual sources. This will lead to better understanding of 
the physical and chemical structure and evolution of disks in the 
early stages of star formation. 
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Appendix A: Snowline test 

Water snowlines were inferred and compared with results from 
Min et al. (2011) in Fig. A. 1 using the minimum mass solar neb¬ 
ula (MMSN) model (2 oc r L5 ). The comparison shows that our 
adopted method reproduces the water snowlines at high accre¬ 
tion rates. For low accretion rates M < 10 M 0 yr 1 , our values 
are slightly smaller yet consistent with those reported in litera¬ 
ture. 



Fig. A.l. Snowlines for the MMSN disk without an envelope: black 
circles show the radii calculated with our method and red squares are 
tabulated values from Min et al. (2011). 


Appendix B: H 2 0, CO 2 , and CO gas fraction 

The H 2 O, CO 2 and CO gas pressure dependent gas fraction 
abundance are shown in Fig. B.l. Figures B.2 and B.3 show the 
midplane C0 2 and CO snowlines as a function of luminosity and 
Ra similar to that of Fig. 9 for H 2 O. 
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Fig. B.l. Gas fraction (n gas /n gas + n ice ) for H 2 0 (top), C0 2 ( middle) and 
CO (bottom) as a function of density and temperatures. 
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Menv = 1 Mg, I/* = 1 I/Q -L* = 5 Lq L * = 15 Lg 



Fig. B.2. Midplane CO 2 snowline as a function of stellar luminosity, accretion rate, and disk radius. The parameters are similar to that of Fig. 9. 


M env — 1 JWq , Z,* — 1 Z/q Z.* — 5 Lg Z,* — 15 Z.q 



Fig. B.3. Midplane CO snowline as a function of stellar luminosity, accretion rate, and disk radius. The parameters are similar to that of Fig. 9. 

























































